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A POLYNOMIAL METHOD FOR DETERMINING LOCAL EMISSION 


INTENSITY BY ABEL INVERSION 
By Chul Park and Dan Moore 
Ames Research Center 


SUMMARY 


The Abel inversion is applied to the transformation of the measured 
line-of-sight integrated radiation intensity from an axially symmetric plasma 
light source into local emission intensity. Consideration of the physical 
features of the plasma light source shows that the most appropriate approxi- 
mation of the radial variation of emission is an even-power polynomial. Thus, 
the measured integrated intensity is represented by a linear combination of 
Abel transforms of even powers. The formula is precise near the axis of sym- 
metry even when the integrated intensity is known at only two or three lateral 
positions. The formula becomes inaccurate at large radial distances and is 
suitable mainly when the emission intensity decreases monotonically toward the 
boundary . 


INTRODUCTION 


The intensity of radiation energy observed from an optically thin plasma 
is an integrated sum of the intensity of emission per unit volume along the 
viewing line. When the plasma is axially symmetric, the radiative emission 
per unit volume, defined here as emission intensity (or simply emission), can 
be calculated from the observed integrated radiation intensity (or integrated 
intensity) by solving Abel's integral equation through inversion (ref. 1). 

In the process of Abel inversion, it is necessary to differentiate the 
experimentally obtained integrated intensity distribution with respect to the 
distance along which it varies. Because the numerical differentiation magni- 
fies the errors in the original data, most methods of Abel inversion result 
in a large magnification of the experimental error, which is proportional to 
the total number of points taken and is severest at the axis of symmetry. 
Typically, there is a tenfold magnification of error when the intensities are 
taken over 20 equally spaced intervals. 

To reduce the effects of experimental error, one ordinarily smooths the 
data by fitting (in a least-squares sense) a polynomial curve through the data 
points. The corresponding curve points are then used as data points before 
applying the Abel inversion to determine the emission. The number of terms 
in the smoothing polynomial can be increased to reduce the sum of squares of 
the deviations of the observed data from the theoretical curve. However, 
high-order polynomials will oscillate, and the variation in slopes will 



increase. Thus, one is compelled to smooth the data with a low-order poly- 
nomial. A data fit with a low-degree polynomial, however, results in an 
analytical error, that is, the error caused by failing to approximate a true 
variation by a simple analytical expression. Therefore, in current Abel inver- 
sion methods, a certain amount of error is introduced regardless of the number 
of points taken to represent the measured distribution of integrated intensity. 

The purpose of the present work is to derive a new method for Abel 
inversion, applicable in certain cases, that will minimize the effects of 
experimental error in determining emission intensity on and near the axis of 
symmetry. Plasma emission intensity is assumed to vary such that it can be 
expressed by a polynomial in even powers of the radial distance. The experi- 
mental data, the integrated intensities, are fitted with a simple function 
that, after Abel inversion, produces a polynomial of even powers. This method 
is applied to two sets of typical experimental data. The proposed method is 
shown to be accurate near the axis of symmetry at least to the same degree as 
existing methods, even when only two or three data points are taken. 


DERIVATION OF METHOD 


X 



Figure 1.- Cross section of axially symmetric, 
optically thin radiation source. 


Figure 1 is a cross section of the 
axisymmetric plasma under consideration. 
When the plasma is optically thin, the 
radiation intensity measured by an 
observer is an integrated sum of volu- 
metric emission that varies with lat- 
eral displacement x of the viewing 
line. It has been shown (ref, 1) that 
integrated intensity P(x) is an Abel 
transform of emission intensity e(r); 
that is. 


P(x) 



e(r)r dr 
\]r 2 - x 2 


( 1 ) 


where r Q is the radius of the plasma boundary. The emission function e(r) 
can be determined (ref. 1) from equation (1) through Abel f s inversion to 
obtain 


e (r) 


Jl f ° (dP/dx) dx 

'J - r 2 

r 


( 2 ) 


In most laboratory plasmas, the following features exist: 

1. At the origin, that is, on the axis of symmetry (r = 0), the emission 
intensity is finite. 
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2. The derivative de/dr approaches zero as r approaches zero. 

3. At the boundary (r = r ) , the derivative (de/dr) r is finite. 

o A 0 

Analytical error can be minimized by fitting the most appropriate 
function to the observed integrated intensity P (x) . Table 1 lists four 
classes of functions that involve polynomials and their Abel inversions. As 
shown in table 1, if an even-power function represents the integrated inten- 
sity 3 the slope of emission function at the boundary is infinite. This fea- 
ture is in conflict with the physical nature of plasma, (3) above, and so the 
method becomes inaccurate near the boundary. The odd-power functions have 
singularities not only at the boundary but also at the origin. For example, 
the linear function P(x) - x implies infinite emission at the origin. For 
other odd-power functions, there is a term of the form r 2n In r in the 
expressions for the emission, which has an infinite curvature at the origin; 
that is, the function e(r) exhibits a finite slope in the region near (but 
not at) the origin, contrary to (2) above. 

Because there is no singularity in the emission function in the range 
0 £ r ^ r Q , the emission can always be approximated by a series in powers of 
r. Table 1 shows the functions P(x) corresponding to the power functions of 
r. For the emission function to be an odd-power function of r, the inte- 
grated intensity must have a logarithmic variation near the origin. Because 
the logarithmic term is undesirable, the preferred form of emission is an 
even-power polynomial 


m 

e ( r ) = ^ b.r 21 
i=0 


(3) 


in which m determines the last term of the series. It can be shown from 
table 1 that the corresponding function for the integrated intensity is 

m 

P(x) = j r^ 2 - x 2 y] a^x 2 ^ (4) 

i=0 


Thus, when the experimentally determined integrated-intensity distribution is 
represented with a function of the form equation (4) , a high degree of accu- 
racy is attained near the origin r = 0. 

The coefficients a^ and b^ in equations (3) and (4) are related through 
Abel f s integral equation (1) (or eq. (2)). This relationship between the a^ 
and b-j_ is derived by substituting equation (3) into equation (1), as 
described in appendix A. For values of m up to 5, explicit formulas for 
determining b^ from given values of a^ are presented in table 2. 


3 



TABLE 1.- ABEL TRANSFORMS OF FUNCTIONS INVOLVING POLYNOMIALS 





^olyfr 2 ) represents an even-power polynomial in r. 
Poly(x 2 ) represents an even-power polynomial in x. 




TABLE 2 - COEFFICIENTS b n IN e(r) =2J b r 2n AS DETERMINED 

n=o n 

m 

FROM THE COEFFICIENTS a n IN P(x) = Jr 2 - X) a n x 2n 


FOR ARBITRARY m 


2 Terms 


= a Q - 0.5air o 2 

= 1.5a! 


3 Terms 


a Q - 0.5air Q 2 - 0.125a2r o 4 
1.5ai - 0.75a2r Q 2 
2 . 125a 2 


4 Terms 

a Q - 0.5ajr o 2 - 0.125a2r o 4 - 0.0625a3r o 6 


1.5ai - 0.75a 2 r o 2 - 0 
2.125a 2 - 0.9375a 3 r o 2 
2.1875ao 


0 . 1875a 3 r 4 


5 Terms 

a Q - 0.5air o 2 - 0.125a2r o 4 - 0.0625a 3 r o 6 - 0.03125a4r 
1.5a! - 0.75a 2 r o 2 - 0.1875a 3 r o 4 - 0.09375a 4 r o 6 
2 . 125a 2 - 0 . 9375a 3 r 0 2 - 0.23437a 4 r o 4 
2 . 1875a 3 - 1 . 09375a 4 r o 2 
2 . 4609a 4 



DETERMINATION OF THE COEFFICIENTS a i 


When the integrated intensity P(xj) known at M + 1 equally spaced 
points x Q = 0, Xi = r 0 /M, x 2 = 2r Q /M, . . . , xj = jr 0 /M, . . . , where M 
is the number of zone divisions, the coefficients a^ of equation (4) can be 
determined from the set of M + 1 linear equations 


j = 0, 1, 2, . . . , M (5) 


However, it must be noted that when j = M, the left side of equation (5) is 
undefined. Physically, this difficulty is caused by the fact that both the 
light-emitting volume and the integrated intensity vanish at x = r Q . By 
taking the limit x + r Q in equation (2), there follows 


PCXj) 


J 1 - Cj/M) 2 


M /. \ 2i 

Z-fe) ■ 


1 = 0 


e Cr 0 ) 


x^r Q Jr Q 2 - : 


( 6 ) 


Because both the numerator and the denominator in the right side of equation 
(6) vanish, it is impossible to determine e(r Q ) experimentally without making 
a suitable assumption on the behavior of P(x) near r Q . Two methods were 
developed to resolve this difficulty to satisfy two types of edge conditions. 


Method 1 

In the first method, the M + 1 equation is deleted so that only M 
linear equations remain to be solved: 


PC*j) 

r 0 <7 1 - (j/M) 2 



i=0 


j = 0, 1, 2, . . . , M - 1 


(7) 

i 


This gives rise to an M term approximation for the emission function. When 
the last equation in (5) is deleted, the emission function is left to take an 
arbitrary finite value at the boundary as dictated by the function values at 
all other points. The e(r Q ) value given by this method may not be correct, 
but it is valid as long as e(r) obeys the polynomial variation law of (3) 
near r 0 . This method should be used if the variation of e(r) with r is 
mild and if e(r 0 ) is known to be finite, as in the case of the radiating 
shock layer around the nose of a body of revolution moving at high speed. The 
solutions for a^ employing this method are presented in table 3. 
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TABLE 3.- COEFFICIENTS b ± IN e(r) 
(xj = jr 0 /M, j = 0, 1 , . . . , M) ; 


m 

= £ b^ 21 AS DETERMINED FROM P(xj) , 
i=0 

METHOD 1 (m = M-l,e(r Q ) t 0 ASSUMED) 


i 


t, 


* 


2 Zones, 2 Terms (M = 2, m = 1) 

b 0 = [1.5P(o) - 1.1547P(r 0 /2)]/r 0 
bi = [-3P (o) + 3.4641P(r 0 /2) ]/r 0 3 

3 Zones, 3 Terms (M = 3, m = 2) 
b 0 = [2 . 0469P (o) - 1 . 3921P (r 0 /3) - 0 . 3144P (2r 0 /3) ] /r 0 

b 1 = [-16 . 0313P (o) + 20 . 2851P (r 0 /3) - 4 . 1507P (2r 0 /3) ] /r Q 3 
b 2 = [18.9844P(o) - 26 . 8480P (r 0 /3) + 8 . 4901P (2r Q /3) ] /r Q 5 

4 Zones, 4 Terms (M = 4, m = 3) 

b 0 = [3 . 2778P (o) - 2 . 7S41P (r 0 /4) + 0 .0770P(r o /2) - 0 . 2688P(3r 0 /4) ]/r Q 
b! = [-43.0000P (o) + 55 . 7710P (r 0 / 4) - 13 . 1636P (r 0 /2) + 0 . 6047P (3r 0 /4) ]/r 0 3 
b 2 = [146 . 6667P (o) - 216 . 8871P (r 0 /4) + 83 . 1384P (r 0 /2) - 13. 1028P(3r o /4) ]/r 0 5 
b 3 = [-124.4444P(o) + 192.7885P(r 0 /4) - 86 . 21 76P (r Q /2) + 1 8 . 8142P (3r Q /4) ] /r Q 7 

5 Zones, 5 Terms (M = 5, m = 4) 

b G = [ 3 . 0694P (o) - 1.6997P(r 0 /5) - 0 . 9633P (2r Q /5) + 0 . 1274P (3r 0 /5) 

- 0.2052P(4r o /5)]/r o 

b i = [ -93 . 2 719P (o) + 124 . 2254P (r 0 /5) - 34 . 1480P (2r 0 /5) + 4 . 1847P(3r 0 /5) 

- 0.8236P(4r o /5)]/r o 3 

b 2 = [579 . 7068P (o) - 877 . 71 89P (r 0 /5) + 369 . 2683P (2r 0 /5) - 78.0015P(3r o /5) 

+ 7.0663P(4r o /S)]/r o 5 

b 3 = [-1260 .9694P(o) + 2010 . 7037P (r Q /5) - 997 . 0724P (2r Q /5) + 283 .9830P (3r 0 /5) 

- 37.4405P(4r o /5)]/r o 7 

b 4 = [834 . 465 IP (o) - 1362 ,6757P(r 0 /5) + 728 . 3808P (2r 0 /5) - 238.4186PC3r D /5) 
j + 39.7364P(4r 0 /5) ]/r D 9 


t 
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AS DETERMINED FROM P(xj) 


m 

TABLE 4.- COEFFICIENTS Ik IN e(r) = £ bir 21 

1 i=0 

(xj = jr 0 /M, j = 0, 1, . . . , M); METHOD 2 (m = M, e(r Q ) = 0 ASSUMED) 

2 Zones, 3 Terms (M = 2, m = 2) 
b Q = [l.SP(o) - 1.1570P(r o /2)]/r o 

b i = [-5 . 25P (o) + 6.9282P(r 0 /2) ]/r 0 3 
b 2 = [3. 75P(o) - 5. 7735P(r 0 /2) ]/r 0 5 

3 Zones, 4 Terms (M = 3, m = 3) 
b 0 = [2 . 2266P (o) "- 1 . 6780P (r 0 / 3) - 0. 1698P(2r 0 /3) ]/r 0 

b x = [-19 . 1016P (o) + 25.1700P(r o /3) - 6.6223P(2r Q /3) ]/r Q 3 
b 2 = [39 . 0234P (o) - 58 . 7300P (r 0 /3) + 24 . 621 2P (2r Q /3) ]/r 0 5 
b 3 = [-22 . 1484P (o) + 35 . 2379P(r 0 /3) - 17. 8291P(2r 0 /3) ]/r 0 7 

4 Zones , 5 Terms (M = 4 , m = 4) 

b 0 = [3.0556P(o) - 2 . 3870P (r 0 / 4) - 0 . 1285P (r 0 /2) - 0 . 1920P (3r 0 /4) ]/r 0 
b x = [-47.9167P(o) + 63.8956P(r 0 /4) - 17. 7054P(r Q /2) + 2 . 3038P (3r Q /4) ]/r Q 3 
b 2 = [200 . 4167P (o) - 305 . 7075P (r 0 /4) + 132. 7906P(r o /2) - 31 .6700P (3r Q /4) ]/r 0 5 
b 3 = [-295 . 5556P (o) + 475 . 5450P(r o /4) - 244 . 2833P(r 0 /2) + 77. 9447P(3r 0 /4) ]/r Q 7 
b 4 = [140 . 0000P (o) - 231 . 3962P (r Q /4) + 129 . 3265P (r Q /2) - 48.3795P(3r 0 /4) ]/r Q 9 

5 Zones, 6 Terms (M = 5, m = 5) 

b 0 = [3 . 7293P (o) - 2 . 822 IP (r 0 /5) - 0 . 2776P(2r 0 /5) - 0 . 1672P (3r 0 /5) 

- 0.11 SOP (4r 0 /5)]/r 0 

b x = [-97. 8760P (o) + 132 . 0571P (r 0 /5) - 38.9323P(2r 0 /5) + 6 . 2401P (3r Q /5) 

- 1 . 4326P (4r 0 /5) ] /r 0 3 

b 2 = [696 . 2967P (o) - 1076 . 0422P (r 0 /5) + 490 . 4206P (2r 0 /5) - 130 .0505P(3r o /5) 

+ 22.4882P(4r 0 /5) ]/r 0 5 

b 3 = [-1937.2940P(o) + 3161 . 1553P (r 0 /5) - 1699 . 8635P (2r 0 /5) + 585 . 9137P (3r 0 /5) 

- 126.9015P(4r o /5)]/r o 7 

b 4 = [2253.0556P(o) - 3775 . 7473P (r 0 /5) + 2202 . 4849P (2r 0 /5) - 871 . 7180P(3r o /5) 

+ 227.3807P(4r o /5)]/r o 9 

b 5 = [-917.9116P(o) + 1561 . 3993P (r 0 /5) - 953 . 8320P (2r 0 /5) + 409 . 7820P (3r 0 /5) 

- 121 . 4169P (4r 0 /5) ]/r 0 1 1 
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Method 2 


In the second method, emission is assumed to be zero at the boundary. 
Thus the full set of M + 1 linear equations in equation (5) was solved with 
the condition 


f 


P( V 

r 0 . h - o/m;i 2 


j = M 


C8) 


This method of solution provides an emission function of M + 1 terms and is 
valid only when the emission decreases at a finite rate to zero at the bound- 
ary (see eq. (6)). The results are presented in table 4. 

Since two methods of solutions of equation (5) were employed, the upper 
limit in the summation in equation (5) will be designated by m, where m = M 
or m = M - 1 depending on which method is used. Appendix B presents the 
method of computing the coefficients a^ , and hence b^, from the given P (x j ) 
for arbitrary values of m up to 5. 


COMPARISON OF RESULTS 


1 


4 


Formulas derived in the preceding 
sections are applied to the two sets of 
experimental data shown in figure 2. 
Data set no. 1 is taken from figure 6 
of reference 2. Data set no. 2 was 
obtained by one of the present authors 
during an experiment with an arc-heated 
wind tunnel (ref. 3) . The measured 
integrated intensity corresponding to 
data set no. 2 is the continuum emis- 

o 

sion at wavelengths near 4000 A from 
ionized nitrogen flowing through a wind 
tunnel. As shown, both sets of data 
have been smoothed so that no error 
would be introduced or amplified 
because of experimental scatter. The 
comparison thus tests the capability of 
the present method to yield the same 
inversion as can be attained by exist- 
ing inversion methods in the absence of 
experimental scatter. Except for the method of reference 2, which uses an 
analog computer, the existing methods use high-order (20- or 30-zone) formulas 
for greater accuracy. The present method limits the number of zones to 5; a 
larger number results in oscillations at larger values of r. 



Figure 2.- Two sets of integrated intensity data 
used for testing the present method. 
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Data Set No. 1 


In figure 3(a), the present 
method of inversion is compared with 
other methods for the first set of data. 
A 5- zone inversion formula of method 1 
(formula 4 of table 3) is chosen for 
this comparison. Corresponding results 
with 4-zone, 4-term and 3-zone, 3-term 
formulas also are shown for the origin 
r = 0. The results of inversions by 
Sirola and Anderson (ref. 2) and 
Hormann (ref. 4) are reproduced from 
reference 2. (The method of Sirola and 
Anderson is based on the analysis of 
Pearce (ref. 5).) Bockasten f s method 
(ref. 6) also was used; the result 
agreed closely (a 2-percent discrepancy) 
with that obtained by Nestor and Olsen 
(ref. 7) . 1 Edels, Hearne, and Young 
(ref. 8) reproduced independently the 
results obtained by Nestor and Olsen 
and by Bockasten. 


Figure 3 indicates that the 
methods of both Hormann and Pearce 
(i.e., Sirola and Anderson) fail at the 
origin. Hormann 1 s method fails because 
of the necessity to evaluate the area 
under an undefined curve near the ori- 
gin. Pearce’s method breaks down because it assumes a step-function variation 
of integrated intensity that cannot describe the emission intensity near the 
origin. The method of Nagler (ref. 9) is similar to that of Pearce. The 
method of Dooley and McGregor (ref. 10) is similar to that of Hormann and has 
the same mathematical difficulty caused by the singularity at the origin. In 
the case of the Hormann inversion, and therefore Dooley and McGregor inversion, 
the error due to the singularity at the origin adds to the deduced value of 
emission over the entire range. 

The methods of Nestor and Olsen and of Bockasten (or equivalently Edels 
et al.) result in definite values for the emission at the origin because the 
linear term is excluded in the integrated intensity (see table 1). However, 
in these methods, the cubic term x^ is allowed as well as x 2 , leading to 
infinite curvature at the origin (see table 1). As a result, the derived 
emission-intensity gradient is finite near the origin as shown, even though 
the gradient of integrated intensity should be zero at the origin. This con- 
tradicts the physical nature of the problem (feature (2) stated earlier) . 



r = r/r 0 


(a) Data set no. 1. 

Figure 3.- Emission intensity obtained by 
inversion of test data in figure 2. 


l The profile inverted by the Nestor-Olsen method shown in figure 3(a) is 

different from the corresponding profile in reference 2. As implied in refer- 
ence 2, the Nestor-Olsen profile was found to be in error. 
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Nevertheless, of the existing methods, these of Nestor and Olsen and of 
Bockasten appear the most accurate; Bockasten f s is preferred because of its 
closer piecewise curve fit (parabolic vs. linear for Nestor and Olsen) to the 
integrated intensity. 

The present method provides a well-defined finite emission value at the 
origin and a smooth parabolic change near the origin. The peak value of 
emission intensity calculated by the present method agrees closely with those 
obtained by the better of the existing methods. The present method fails, 
however, at large values of r* = r/r Q , that is, r" > 0.6 for the present test 
data, because: 


1. The formula used here (method 1) leaves the emission arbitrary at 
the boundary. 

2. Expansions of the form of equation (4) tend to oscillate as r/r 0 1. 

Point (2) is a result of the emission-function peak at an intermediate point, 
around r = 0.5, which decays rapidly as "r increases toward unity. Thus, 
the present method is not suitable for calculating the emission at large radial 
distances when the slopes are changing rapidly (see fig. 2). At the origin, 
however, even a 3- zone formula of either method 1 or 2 gives a satisfactory 
result. As seen from figure 3, existing methods are satisfactory at large __ 
radial distances. It is possible to improve the present method for large r, 
for example, by expanding in power series around r Q instead of the origin; 



(t>) Data set no. 2. 
Figure 3-- Concluded. 


however, no such attempt was made since 
most existing methods yield satisfactory 
results in this region. 

Data Set No. 2 

Figure 3(b) compares different 
formulas of methods 1 and 2 with that of 
Bockasten for the second set of data. 

To avoid congestion, results are shown 
for selected radii for 2-, 4-, and 
5-zone cases. For the 2-zone cases, the 
emission values are shown only at the 
origin and approximate midradius 
Y = 0.47, because these are sufficient 
to define the curve. 

The Bockasten inversion, like that 
of Nestor and Olsen (fig. 3(a)), yields 
a sharp peak in emission at the origin. 
For the reason noted earlier, the 
Bockasten value at the origin cannot be 
trusted. For the rest of the range, 
except near the plasma boundary, 
methods 1 and 2 formulas of order 3 and 
above yield results that agree closely 
with that of Bockasten. Generally, the 
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higher order formulas tend to oscillate near the boundary. The amplitude of 
the oscillations are greater for method 1 than for method 2, which assumes 
that emission vanishes at the boundary. In particular, the 3-zone formula of 
method 2 agrees closely with the Bockasten result over the entire radius 
range. This agreement may be fortuitous and cannot be expected as a rule. 
Because none of the existing methods is accurate at the origin, it is not 
possible to assess the absolute accuracy of the present method at the origin. 
Such accuracy is strongly suggested, however, by the consistency of results 
yielded by all the formulas of different order, and by the mathematical rigor 
of the present method. Even a 2-zone, 2-term formula gives a value accurate 
to within 5 percent at the origin. 

The present method results in high accuracy for data set no. 2 because 
the emission variation with radius is well behaved (fig- 3(b)). That is, the 
emission decreases monotonically toward the boundary, and the change in the 
slope of decay is gradual over the entire radius range. Because such a varia- 
tion can be approximated accurately by a polynomial, a high degree of accuracy 
can be expected. 


DISCUSSION 


The sensitivity of the inversion to the initial input error in integrated 
intensity is shown in table 5. The numbers in the table represent the ratio 
of the output error in emission at the radii indicated to the input error in 
integrated intensity. Thus, for example, if one uses the 4-zone, 5-term for- 
mula for data set no. 1, and if the integrated intensity at r = 0.5r o is 
misread by 1 percent, the emission suffers an error of 1.18 percent at that 
point. The 2-zone, 2-term, and 3-term calculation for data set no. 1 are not 
listed because the 2- term formulas are inappropriate. 

As shown, the present method is much less sensitive to experimental error 
than the existing methods, at least in the range 0 < r < 0.5r o . Note that 
the amplification of initial error is small with the present method even though 
the coefficients in the expressions of b^ in tables 3 and 4 are large num- 
bers, especially for large m. It is essential, however, that the values of 
P(jr 0 /M) be represented by a sufficiently large number of digits in the calcu- 
lation of b^ for all i. A lower order formula magnifies the input error 
to a lesser degree. Again, this is expected: the greater the number of inter- 

vals at which a value of integrated intensity is read, the smaller the magnifi- 
cation of error by differentiation. For the same zone division, method 1 
magnifies the error the least. 

Any error in the determination of the overall radius of r Q is magnified 
by the inversion process. To determine the factor of magnification of such an 
error, one first changes r 0 deliberately by an amount Ar 0 so that the jth 
zone point is shifted by (j/M)Ar 0 , and computing the emission function e on 
the basis of the values of P at the new zone points. This calculation is 
carried out for data set no. 2 for the present method and for the 10-zone 
Bockasten formula. Table 6 shows the results of the calculation. The figures 
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TABLE 5.- COMPARISON OF INVERSION METHODS FOR MAGNIFICATION 3 

OF EXPERIMENTAL ERROR 



Data set no. 1 

Data 

set no. 2 


r = 0 

r = r 0 /2 

r = 0 

r = r 0 /2 

Bock as ten (ref. 6) 
M = 20 

18.8 

1.88 

14.2 

3.08 

Nestor and Olsen 
(ref. 7) 





M = 20 

15.8 

2.17 

11.8 

3.40 

M = 50 

39.5 

5.43 

29.5 

8.50 

Method 1 





m = 2 

— 

— 

1.44 

.20 

m = 3 

2.77 

.61 

1.90 

.27 

m = 4 

3.47 

.76 

2.38 

.33 

m = 5 

4.16 

.91 

2.85 

.40 

Method 2 





m = 2 

— 

— 

1.42 

.26 

m = 3 

3.11 

.91 

2.03 

.37 

m = 4 

4.02 

1.18 

2.85 

.52 

m = 5 

4.97 

1.46 

3.46 

.63 


cl 

Ratio of the error in emission intensity e(r) to the error in the 
integrated intensity P(x) . 


TABLE 6.- COMPARISON OF PRESENT AND BOCKASTEN INVERSION METHODS 


FOR MAGNIFICATION OF ERROR IN r Q * 



M = 2 

M = 3 

M = 4 

M = 5 

M = 10 

Method 1 

1.14 

1.51 

1.67 

1.72 

— 

Method 2 

1.14 

1.67 

1.70 

1.57 

— 

Bock as ten 

— 

— 

— 

— 

1.18 


*Ratio of error in emission intensity at the origin caused by the error 
in r Q to the error in r 0 , for data set no. 2. 
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in the table represent the magnification factor, that is, the ratio of the 
error in the resulting emission intensity at the origin e(0) caused by the 
error in r Q to the error in r Q itself. The 10-zone (rather than 20-zone) 
Bockasten formula was used because it was difficult to read small changes in 
P(x) caused by the shifting of x by (j/M)Ar 0 . As indicated by table 6, the 
magnification of the error in r G by the inversion remains the order of unity 
for both the present and Bockasten T s methods. 

In the present method, a truncation error exists in addition to the 
experiment-originated errors mentioned above. Because the basic polynomial 
representations, equations (3) and (4), truncate a Maclaurin infinite series 
at the mth term, there always exists an error of the same magnitude as the 
m + 1 term. The accuracy of the present method depends, among other things, 
on how accurately the functions e and P can be approximated by the forms of 
equations (3) and (4). The truncation error is small near the origin, but 
becomes large at large r, especially for large m, making the present method 
invalid at large r. 


The curvature of emission variation with radius at the axis is given in 
the present method as 


d 2 e 

dr 2 


2b i 


where b\ is as listed in tables 3 and 4. The accuracy of the derived curva- 
ture will improve as the degree of the polynomial is increased provided that 
the input data are absolutely accurate. None of the existing methods, however, 
can satisfactorily resolve emission profile curvature on the axis of symmetry. 

When the integrated intensity is not known at equally spaced points 
x = jr 0 /M, the formulas of tables 3 and 4 cannot be used. The present method 
can still be applied, however, with a minor modification: The curve of the 

form of equation (4) must be fitted to the experimental data; then the coeffi- 
cients b^ can be determined as in appendix A or table 2. 

Method 2 should be used if the emission is known to vanish at the 
boundary and if the values near the boundary are desired. Otherwise, method 1 
is preferable because it is insensitive to initial input error. The required 
number of zones depends on the type of variation of integrated intensity. If 
the integrated intensity varies such that the emission varies monotonically 
with radius, and if there are no more than two inflection points, as in fig- 
ure 4, a 3-zone subdivision should be adequate. If there is no inflection 
point in integrated intensity, as when the boundary is a bow shock wave over a 
spherical body moving through the atmosphere at very high speeds, a 2- zone 
formula may be accurate enough. The order of zoning should be increased as 
the number of inflection points is increased, because otherwise the polynomial 
cannot properly describe the variation. If it is not known a priori how many 
inflection points the emission function has, the calculation should be 
repeated with various zonings. The required number of zones will have been 
reached when further increases fail to change the calculated emission values. 
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I 


The present method can be used most advantageously in experiments to 
determine the conditions of plasma produced in facilities for aerothermody- 
namic testing. In most such experiments in wind tunnels or shock tubes, it 
is necessary only to know the condition of the plasma at and near the center- 
line, that is, the emission and curvature at the origin. Also, in many cases, 
the flows produced in such facilities have maximum emission on the axis of 
symmetry, which is favorable for the present method. 


CONCLUSIONS 


The method for Abel inversion described uses only a few measured values 
of integrated radiation intensity. The emission intensity is described by an 
even-power polynomial in the radial coordinate r, the coefficients of which 
are determined to best fit the measured intensity. The inversion is accurate 
near the axis of symmetry, and breaks down at large radial distances. It is 
also more accurate when emissivity decreases monotonically toward the boundary 
than otherwise. The magnification of error in the inversion is much smaller 
than that present in other methods. 

Ames Research Center 

National Aeronautics and Space Administration 
Moffett Field, Calif., 94035, Oct. 2, 1969 


15 


APPENDIX A 


DETERMINATION OF b j_ FROM ai 


By differentiating equation (4) , one obtains 

m 


al = T7—^ S x2ktl [ 2(k * 15 Vi r o 2 - ( 2k * » vl 


k=0 


Then one evaluates the integral 




f 

Jr 


x 2k+1 dx 




x 2 - r 2 ) (r 2 - x 2 


_ _ I JL \ ^ /2«,\/2k - 2 A zt 2(k- 
7" 2 4 k W\ k - £ / ° 


JO 


£=0 


Substitution of equations (Al) and (A2) into equation (2) yields 


m 


= OZ [ (2k * « Vi r o 2 - t2k + 1Ja k] Q i. 


k=o 


(Al) 


(A2) 


(A3) 


Equating the coefficients of corresponding powers of r in expressions (3) 
and (A3) , we obtain the coefficients 


m 


b. = 
1 


^2 2 ^ ( 2 k - + 1)a k - C2k - 2)a k + l r o 2 ] 

k=i 


(A4) 


In terms of the a^, equation (A4) becomes 


b i [ (2k + i)B k,i ■ 2 kB k-i,ii r o 2k - 2 s 


where 


B 


k,i 


1 j_ / 2k ■ 

2 4 k V k - 



for k ^ i 


(A5) 


(A6) 


The explicit numerical form of equation (A5) is presented in table 2 for values 
of m up to 5. 
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APPENDIX B 


DETERMINATION OF a ± AND bi FROM P(xj) FOR EQUALLY SPACED x- 


In the solution of equation (5), the values for P(x) are assumed known 
at equally spaced points xj 

X 3 = M r o > j = 0, 1, 2, . . . , M 

Then equation (5) can be written as 


P(jr Q /M) 


r Jl - (j/M) 2 


m \2k 

k=o ' 


(Bl) 


Letting Y be a vector with components yj , X a matrix with components 
Xj^j and A a vector with components a^, where 


y • = 


P(jr Q /M) 


- 1 r 1 - (j/M) 2 


X. v = 


m • 


and X =1 


we have the matrix representation of equation (Bl) 

Y = XA 

Now we can find a solution matrix X" 1 so that the a k can be given by 

m 

“k ■ E x kjq 

j=o 


(B2) 


Substitution of equation (B2) into (A5) yields 
m m 

»i ‘ E E [ (2k * 1)B k,i - 2kB k-l,i] r o 2k_2nX kj y ) <B3) 

k=i j =0 

where Bj. ^ is given in equation (A6) . The numerical form of equation (B3) 
is presented in tables 3 and 4 for values of m up to 5. 
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